In [337]:
def embed_vid(outvid):
    video = io.open(outvid, 'r+b').read()
    encoded = base64.b64encode(video)
    return HTML(data='''<video alt="test" width="950" height="500" loop="true" controls>
                <source src="data:video/mp4;base64,{0}" type="video/mp4" loop="true" />
             </video>'''.format(encoded.decode('ascii')))

Storm Track Analysis

The tracking algorithm is a fork of the Tint (Tint is not Titan) tracking algorithm (http://openradarscience.org/TINT/). The framework has been modified that it can be also applied to model output data that is not not stored in Py-ART radar data container.

  • The analysis is based averages of one week worth of strome cell tracking:

1. Median of strom area, duration, avg., max. rain-rates and # of storm cells

In [17]:
#Plot the Medians
um044_n = len(UM044_ens.coords['dim_0'])
um133_n = len(UM133_ens.coords['dim_0'])
cpol_n = len(OBS.dataset['total'].coords['dim_0'])
cpol_n2 = len(CPOL.coords['dim_0'])
CPOL_pdf2 = ravel(OBS.dataset['obs'])
var=['avg_area','dur', 'avg_mean', 'max_mean', 'dist', 'v']
names=['area', 'duration', 'avg-rain', 'max-rain', 'distance', 'speed', '# storms']
#print(CPOL_pdf[var].median())
medians = pd.DataFrame({'a': list(CPOL_pdf[var].median().values)+[int(cpol_n)],
                        'b': list(CPOL_pdf2[var].median().values)+[int(cpol_n2)],
                        'd': list(UM044_pdf[var].median().values)+[int(um044_n)],
                        'c': list(UM133_pdf[var].median().values)+[int(um133_n)]} )
                       #index=('Area','Duration', 'Mean-Rain', 'Max-Rain'))
medians.index=names
medians.columns=['Bootstrap', 'Cpol', 'UM 1.33km', 'UM 0.44km']
#print('Medians:')
                 
medians.round(2)
Out[17]:
Bootstrap Cpol UM 1.33km UM 0.44km
area 107.66 119.38 75.52 57.92
duration 60.00 74.67 60.00 50.00
avg-rain 4.52 4.50 4.78 5.55
max-rain 6.75 6.60 6.88 8.26
distance 14.80 16.42 10.02 10.62
speed 12.67 12.83 10.03 12.78
# storms 1184.00 42.00 50.00 73.00

2. 2D-Histograms (storm area vs duration):

In [18]:
#Create Hex-Bin plot
mpld3.disable_notebook()
var=['avg_area','dur', 'max_mean', 'avg_mean']
medians = pd.DataFrame({ 'Bootstrap':CPOL_pdf[var].median(),
                        'CPOL':CPOL_pdf2[var].median(),
                        'UM 1.33km':UM133_pdf[var].median(),
                        'UM 0.44km':UM044_pdf[var].median()})
#medians.loc['avg_area'] /= 2.5**2
histdata =  [UM044_pdf[var].dropna(), UM133_pdf[var].dropna(), CPOL_pdf2[var].dropna(), CPOL_pdf[var]][::-1]
titles = ['UM 0.44km ens', 'UM 1.33km ens', 'CPOL', 'Bootstrap'][::-1]
fig = plt.figure(figsize=(10,8))
colm = colmap2.Blues
colm.set_under('w', alpha=0)
fig, (ax1, ax2, ax3, ax4) = plt.subplots(1, 4, sharey=True)
fig.subplots_adjust(bottom=0.07, right=0.98, left=0.01, top=0.35, wspace=0.01)
cbar_ax = fig.add_axes([0.15, 0.02, 0.7, 0.01])
hexbin_data = []
gridsize = 10
vmin=0.0001
vmax=0.1
nticks=10
YMax, XMax = 260, 70*2.5**2
#YMax, XMax = 50, 70*2.5**2
for i, ax in enumerate((ax1, ax2, ax3, ax3)):
    data = histdata[i][var[:2]]
    #data[var[0]] /= 2.5**2
    #data = data.loc[(data[var[0]] <= XMax) & (data[var[1]] <=YMax)]
    X = data[var[0]].values
    Y = data[var[1]].values
    ax.set_ylim(0,YMax)
    ax.set_xlim(0,XMax)
    hb =  ax.hexbin(X, Y,  gridsize=gridsize, extent=(0,XMax,0,YMax))
    hexbin_data.append(np.ones_like(Y, dtype=np.float) / hb.get_array().sum())
    
plt.cla()
medians = OrderedDict()

for i, ax in enumerate((ax1, ax2, ax3, ax4)):
    data = histdata[i][var[:2]]
    #data[var[0]] /= 2.5**2
    #data = data.loc[(data[var[0]] <= XMax) & (data[var[1]] <=YMax)]
    X = data[var[0]].values
    Y = data[var[1]].values
    ax.set_ylim(0,YMax)
    ax.set_xlim(0,XMax)
    im = ax.hexbin(X, Y,  gridsize=gridsize, cmap=colm, marginals=False, extent=(0,XMax,0,YMax),
                     vmin=vmin, vmax=vmax, C=hexbin_data[i], reduce_C_function=np.sum)
    
    ax.set_title(titles[i], fontsize=24)
    ax.grid(color='w', linestyle='', linewidth=0)
    ax.tick_params(labelsize=24)
    ax.xaxis.set_ticks(ax.xaxis.get_ticklocs()[:-1])
    x, y = histdata[i][var[0]].median(), histdata[i][var[1]].median()
    z, zz=  histdata[i][var[2]].median(), histdata[i][var[-1]].median()
    sx, sy = histdata[i][var[0]].std(), histdata[i][var[1]].std()
    medians[titles[i]] = '%2.1f km^2, %2i min (max: %2.1f mm/h, mean: %2.1f mm/h);'%(x,y, z, zz)
    ax.hlines(y*np.ones_like(histdata[i][var[0]]),0,x, 'firebrick', lw=3)
    ax.vlines(x*np.ones_like(histdata[i][var[1]]),0,y, 'firebrick', lw=3)
    ax.scatter([x], [y], marker='o', s=400, c='firebrick', alpha=0.8)
    if i == 0:
        ax.set_ylabel('Duration [min]', fontsize=24)
    ax.set_xlabel('Rainfall Area [km$^2$]', fontsize=24)
ary=im.get_array()/im.get_array().sum() * 100
im.set_array = ary
cbar = fig.colorbar(im, cax=cbar_ax, orientation='horizontal')
cbar.ax.tick_params(labelsize=24)
cbar.set_label('Density [ ]',size=24)
#cbar.set_ticks(np.linspace(vmin, vmax, nticks).round(2))
#cbar.set_ticklabels(np.linspace(vmin, vmax, nticks).round(2))
#print('Medians:')
#medians
#print(' '.join(['%s: %s'%(k, v) for (k,v) in medians.items()]))
plt.show()
<matplotlib.figure.Figure at 0x7f3b7c4b8b70>
In [19]:
colm = matplotlib_to_plotly(CubeYF_20.get_mpl_colormap(N=8, gamma=2.0),8, rgb=False)

3. Storm tracks above the 80th percentiles

  • marks start of track, colors for different ensemble member
In [21]:
#Plot the tracks
mpld3.enable_notebook()
fig = plt.figure(figsize=(10,7))
fig.subplots_adjust(right=0.94, bottom=0.45, top=0.95,left=0.01, hspace=0, wspace=0)
#cbar_ax = fig.add_axes([0.12, 0.17, 0.74, 0.02])
o = namedtuple('Sim', 'dataset percentiles')
tmp_obs = o({'obs': CPOL_t.dataset['obs']}, {'obs': CPOL_t.percentiles['obs']})

with nc(CPOLF) as fnc:
    lon=fnc.variables['longitude'][:]
    lat=fnc.variables['latitude'][:]
tp = 'mean'
num=80
titels = ['CPOL', 'UM 1.33km', 'UM 0.44km']
m = None
for i, data in enumerate(((tmp_obs, storm_obs), (UM133_t, storm_UM133), (UM044_t, storm_UM044))):
    tracks, stormId = data
    ax = fig.add_subplot(2,3,i+1)
    ax2 = fig.add_subplot(2,3,i+4)
    ax.set_title(titels[i], fontsize=18)
    for ii, tr in enumerate(tracks.dataset.keys()):
        if ii == 0:
            draw_map = None
            m2 = None
        else:
            draw_map = m
        Id = [stormId.Sim[tr]]
        perc = tracks.percentiles[tr][tp][num]
        ax, m, im = plot_traj(tracks.dataset[tr], lon, lat, ax=ax, mintrace=2, thresh=('mean', perc), 
                  color=colm[ii], create_map=draw_map, basemap_res='f', lw=0.5, size=20, particles=None)
        ax2, m2, im = plot_traj(tracks.dataset[tr], lon, lat, ax=ax2, mintrace=2, thresh=('mean', perc),
                  color=colm[ii], create_map=m2, basemap_res='f', lw=0.5, size=20, particles=Id)
        #break
plt.show()

4. Classification of strom area, rain-rate and duration by rain-rate quintiles

In [1115]:
mpld3.disable_notebook()
variables=dict(dur=('Duration [min]',300), avg_area=('Area [km$^2$]', 100*2.5**2), avg_mean=('Rain-rate [mm/h]',12))
fig = plt.figure()
fig.subplots_adjust(right=0.94, bottom=0.025, top=0.6,left=0.01, hspace=0, wspace=0)
i = 0
for y, prop in variables.items():
    label, ylim = prop
    npl = len(list(variables.keys()))
    i += 1
    ax = fig.add_subplot(npl,1,i)
    ax = sns.boxplot(x="quant", y=y, hue="run", data=df_stack, palette="muted", ax=ax)
    #ax = sns.stripplot(x="quant", y=y, hue="run", data=df_stack, jitter=True, palette="Set2", dodge=True)
    ax.legend(loc=0, fontsize=24)
    ax.tick_params(labelsize=24)
    ax.set_xlim(0.5,5.5)
    ax.set_ylim(0,ylim)
    ax.yaxis.set_ticks(ax.yaxis.get_ticklocs()[1:])
    ax.set_xlabel('Quintiles', fontsize=26)
    ax.set_ylabel(label, fontsize=26)
    #break
#fig.savefig(os.path.join(os.environ['HOME'], 'Todds_Rainfall_2.pdf'), bbox_inches='tight', format='pdf', dpi=300)

5. Stormtrack Percentiles vs Rain-rate

In [516]:
P = list(np.arange(0,99))+[99, 99.9, 99.99, 99.999, 100]

PERC = pd.DataFrame({'Obs':np.percentile(CPOL_pdf2['avg_mean'].values, P), 
                     'UM 1.33km': np.percentile(UM133_pdf['avg_mean'].values, P),
                      'UM 0.44km': np.percentile(UM044_pdf['avg_mean'].values, P)},index=P)

from mpl_toolkits.axes_grid.inset_locator import inset_axes
#Plot the percentages
fig=plt.figure()
ax = fig.add_subplot(111)
plt.subplots_adjust(bottom=0.05, right=0.9, top=0.4)
ax.plot(PERC.index, PERC['Obs'].values, color='lightseagreen', linestyle='-', label='CPOL',lw=3)
ax.plot(PERC.index, PERC['UM 1.33km'].values, color='fuchsia', linestyle='--', label='UM 1.33km',lw=3)
ax.plot(PERC.index, PERC['UM 0.44km'].values, color='fuchsia', linestyle='-', label='UM 0.44km', lw=3)
ax.fill_between(PERC.index, member.min(axis=0)-0.25, member.max(axis=0)+0.25, color='cornflowerblue', alpha=0.2)
ax.plot(PERC.index, member.mean(axis=0), color='lightseagreen', linestyle='--', label='Bootstrap',lw=3)
ax.set_xlim(1,100)
#ax.set_ylim(10,100)
#ax.set_yscale('log')
#ax.set_xscale('log')
ax.set_xlabel('Percentile [\%]', fontsize=26)
ax.set_ylabel('Rain-rate [mm/h]', fontsize=26)
ax.legend(loc=0, fontsize=24)
ax.tick_params(labelsize=24)
ax.grid(color='r', linestyle='-', linewidth=0)
In [27]:
#Bootstrapping stuff:

6. Individual Cases (12/11 2006)

In [331]:
#Now make a hovemoller plot of the date
mpld3.disable_notebook()
cmap = colmap2.Blues
cmap.set_under('w')
cmap.set_bad('w')
import matplotlib.gridspec as gridspec
from itertools import product
fontsize=16

## Define start and end time
start_time = pd.DatetimeIndex([datetime(2006,11,12,12,0, tzinfo=timezone)])[0]
end_time = pd.DatetimeIndex([datetime(2006,11,12,18,0, tzinfo=timezone)])[0]
hfmt = dates.DateFormatter('%H')
## Get the information for cpol obs
dy1 = len(OBS_grid.dataset[1].coords['x'])//2
obs = OBS_grid.dataset[1].variables['lsrain'][...,:-20,:].mean(dim=('x'))
obs_tsteps = pd.DatetimeIndex(OBS_grid.dataset[1].coords['t'].values).tz_localize(utc).tz_convert(timezone)
obs_lon = OBS_grid.dataset[1].variables['longitude'][0,:]
sIdx = np.fabs((obs_tsteps - start_time).total_seconds().values).argmin()
eIdx = np.fabs((obs_tsteps - end_time).total_seconds().values).argmin()+1
obs_tsteps = obs_tsteps[sIdx:eIdx]
obs_data = np.ma.masked_less_equal(obs[sIdx:eIdx].values, 0.000)
#sIdx, eIdx = 0, -1

## Get the information from the simulations
### 044 SIM
dy1 = len(UM044ens_grid.coords['lat'])//2
um_tsteps1 = pd.DatetimeIndex(UM044ens_grid.coords['t'].values).tz_localize(utc).tz_convert(timezone) 
um1 = UM044ens_grid.variables['lsrain'][...,:-20,:].mean(dim=('lat','surface'))
sIdx = np.fabs((um_tsteps1 - start_time).total_seconds().values).argmin()
eIdx = np.fabs((um_tsteps1 - end_time).total_seconds().values).argmin()+1
um_data1 = np.ma.masked_less_equal(um1[:,sIdx:eIdx].values, 0.00)
um_lon1 = UM044ens_grid.coords['lon'][:]
um_tsteps1 = um_tsteps1[sIdx:eIdx]
tit = ['Obs']
## 133 SIM
dy1 = len(UM133ens_grid.coords['lat'])//2
um_tsteps2 = pd.DatetimeIndex(UM133ens_grid.coords['t'].values).tz_localize(utc).tz_convert(timezone) 
um2 = UM133ens_grid.variables['lsrain'][...,:-20,:].mean(dim=('lat','surface'))
sIdx = np.fabs((um_tsteps2 - start_time).total_seconds().values).argmin()
eIdx = np.fabs((um_tsteps2 - end_time).total_seconds().values).argmin()+1
um_data2 = np.ma.masked_less_equal(um2[:,sIdx:eIdx].values, 0.00)
um_lon2 = UM133ens_grid.coords['lon'][:]
um_tsteps2 = um_tsteps2[sIdx:eIdx]
## Stack the simulations together
um_data = [(um_data2, um_lon2, um_tsteps2) , (um_data1, um_lon1, um_tsteps1)]

#fig, ax = plt.subplots(3,3, sharey=True, sharex=True)
fig = plt.figure(figsize=(15,10))
#ax=ax.ravel()


outer_grid = gridspec.GridSpec(3, 3, wspace=0.0, hspace=0.12)
inner_grid = gridspec.GridSpecFromSubplotSpec(1, 2, subplot_spec=outer_grid[0], wspace=0.0, hspace=0.0)
ax = plt.Subplot(fig, inner_grid[0])
im = ax.pcolormesh(obs_lon, obs_tsteps.tz_localize(None), obs_data,vmin=0.0,vmax=5,cmap=cmap)

fig.add_subplot(ax)

ax.set_title('Obs', fontsize=fontsize)
ax.tick_params(labelsize=fontsize-3)
#ax.xaxis.set_ticks(list(ax.xaxis.get_ticklocs()[:-1][::3]))
ax.xaxis.set_ticks([])
#ax.set_ylabel('Local Time', fontsize=28)
ax.yaxis.set_major_formatter(hfmt)
add = ('UM 1.33km', 'UM 0.44km')
for i in range(1,9):
    inner_grid = gridspec.GridSpecFromSubplotSpec(1, 2, subplot_spec=outer_grid[i], wspace=0.0, hspace=0.0)
    for j, data in enumerate(um_data):

        um, lon, tsteps = data
        ax = plt.Subplot(fig, inner_grid[j])
        m = ax.pcolormesh(lon, tsteps.tz_localize(None), um[i-1],vmin=0.0,vmax=5,cmap=cmap)
        if i in (8,7,6):
            ax.tick_params(labelsize=fontsize-4)
            ticks =[' ']+list(ax.xaxis.get_ticklocs()[1:][::3].round(1).astype(str))
            #ax.xaxis.set_ticks(ticks)
            ax.xaxis.set_ticks(list(ax.xaxis.get_ticklocs()[:-1][::3]))
        else:
            ax.set_xticks([])
        tit = add[j]+ r' (\rom{{{}}})'.format(i)
        #tit = add[j]+datetime.strptime(ensembles[i-1], '%Y%m%dT%H%MZ').strftime(' (%e %b %H UTC)')
        #if j == 0:
            
        if i == 3 or i == 6:
            if j == 0:
                #ax.set_ylabel('Local Time', fontsize=28)
                ax.tick_params(labelsize=fontsize-3)
                ax.yaxis.set_major_formatter(hfmt)
            else:
                ax.set_yticks([])
        else:
            ax.set_yticks([])
        ax.set_title(tit, fontsize=fontsize-4)


        fig.add_subplot(ax)
fig.subplots_adjust(right=0.98, bottom=0.25, top=0.99,left=0.01, hspace=0.1, wspace=0)
cbar_ax = fig.add_axes([0.14, 0.18, 0.74, 0.01])
cbar=fig.colorbar(im, cax=cbar_ax, orientation='horizontal')
cbar.ax.tick_params(labelsize=fontsize-4)
cbar.set_label('Avg. Rain-rate [mm/h]',size=fontsize-2)
fig.text(0.5, 0.195, 'Longitude', ha='center', fontsize=fontsize)
fig.text(-0.03, 1-0.38, 'Local Time (12. Nov. 2006)', va='center', rotation='vertical', fontsize=fontsize)
#for i in (6, 7, 8):
#    ax[i].set_xlabel('Longitude', fontsize=28)
Out[331]:
<matplotlib.text.Text at 0x7f3b1b60a518>
  • Ensemble member III and V:
In [33]:
um044_n = len(UM044_trackData.coords['dim_0'])
um133_n = len(UM133_trackData.coords['dim_0'])
cpol_n = len(CPOL_trackData.coords['dim_0'])

cpol_d = round(CPOL_trackData.to_dataframe()['dist'].max(),2)*2
um044_d = np.array(UM044_trackData['dist'].max(axis=1).median()).round(2)*2
um133_d = np.array(UM133_trackData['dist'].max(axis=1).median()).round(2)*2


var=['avg_area','dur', 'dist', 'avg_mean', 'max_mean']
names=['area', 'duration', 'distance', 'avg-rain', 'max-rain', 'distance','# storms']
#print(CPOL_pdf[var].median())
medians = pd.DataFrame({'b': list(CPOL_trackData[var].to_dataframe().median().values.round(2))+[cpol_d,int(cpol_n)],
                        'd': list(UM044_trackData[var].to_dataframe().median().values.round(2))+[um044_d,int(um044_n)],
                        'c': list(UM133_trackData[var].to_dataframe().median().values.round(2))+[um133_d,int(um133_n)]} )
                       #index=('Area','Duration', 'Mean-Rain', 'Max-Rain'))
medians.index=names
medians.columns=['Cpol', 'UM 1.33km', 'UM 0.44km']
#print('Medians:')               
medians.round(2)
Out[33]:
Cpol UM 1.33km UM 0.44km
area 99.31 83.37 58.33
duration 121.60 70.00 50.00
distance 16.72 11.64 9.13
avg-rain 4.69 5.47 6.21
max-rain 7.09 7.86 10.21
distance 75.76 58.78 46.40
# storms 3.00 5.00 7.00

8. Cold-Pools

Cold-Pools of ensemble member III and V for a storm on 12/11 2006

  • Cold-Pools along the storm track
In [518]:
fig=plt.figure(figsize=(15,13))
fontsize=20
nplot=1
cold_pool_times = {}
hfmt = dates.DateFormatter('%H:%M')
for n, i in enumerate((3, 4)):
    for j, (name, data, coords) in enumerate(zip(('UM133', 'UM044'), (UM133_thetae, UM044_thetae), (UM133_surf, UM044_surf) )):
        pax = fig.add_subplot(2,2,nplot)
        data -= data.mean()
        if i == 4:
            tit = ('UM 1.33km','UM 0.44km')[j]+ r' (\rom{{{}}})'.format(i+1)
        else:
            tit = ('UM 1.33km','UM 0.44km')[j]+ r' (\rom{{{}}})'.format(i)
        time = pd.DatetimeIndex(coords.coords['t'].values[idx[i][name]['time']]).tz_localize(utc).tz_convert(timezone)
        time = time.tz_localize(None)
        tsurf = get_data(data, idx, i, name)[...,0].T
        
        X = 5 * np.arange(tsurf.shape[1])
        diff = np.fabs(np.diff(tsurf, axis=0))
        peak = np.unravel_index(diff.argmax(), diff.shape)
        
        im = pax.pcolormesh(X, time, tsurf, vmin=-1, vmax=10, cmap=colmap.GMT_haxby, shading='gouraud')
        pax.tick_params(labelsize=fontsize-2)
        pax.plot([X[0],X[-1]], [time[0],time[-1]])
        #pax.scatter(X[peak[1]], time[peak[0]+1], c='k')
        cold_pool_times[i] = peak[0]+1
        pax.set_title(tit, fontsize=fontsize)
        pax.yaxis.set_major_formatter(hfmt)
        nplot += 1   
    
fig.subplots_adjust(right=0.98, bottom=0.3, top=0.99,left=0.01, hspace=0.2, wspace=0.15)
cbar_ax = fig.add_axes([0.14, 0.25, 0.74, 0.01])
cbar=fig.colorbar(im, cax=cbar_ax, orientation='horizontal')
cbar.ax.tick_params(labelsize=fontsize-2)
cbar.set_label('Density Potential Temperature Pertubation [$^\circ$C]',size=fontsize)
  • Snapshot:
In [1379]:
fig=plt.figure(figsize=(15,10))
fontsize=14
ds = 10
hfmt = dates.DateFormatter('%H:%M')

#
first = True
M = []
plev = 0
pdata = 'pres'
if pdata == 'surf':
    coldpool_co = {}
for tidx in range(111,112):
    nplot = 1
    sys.stdout.flush()
    #sys.stdout.write('\r Creating Plot %i of 114'%tidx)
    sys.stdout.flush()
    for n, i in enumerate((3, 5)):
        if pdata == 'pres':
            UM133_plev, UM044_plev = get_file('vert_cent', num=i-1, UMdir=os.path.join(UMdir, 'tmp'),
                                              first='20061112_0000', last='20061112_1200', remap_res='native')
            UM133_rain, UM044_rain = get_file('rain', num=i-1, UMdir=os.path.join(UMdir, 'tmp'),
                                              first='20061112_0000', last='20061112_1200', remap_res='native')
            P = UM133_plev['p'].values[plev]
            UM133_thetap = calc_thetap(UM133_plev['temp'][:, plev, :], P,
                                       UM133_plev['q'][:, plev, :].values, UM133_rain['lsrain'][:,0,:].values)
            UM044_thetap = calc_thetap(UM044_plev['temp'][:, plev, :], P,
                                       UM044_plev['q'][:, plev, :].values, UM044_rain['lsrain'][:,0,:].values)
            #UM133_plev, UM044_plev = get_file('vert_cent', num=n)
            #UM133_thetap = calc_thetap(UM133_plev['temp'][:, plev, :].values, UM133_plev['p'][plev].values,
            #                           UM133_plev['q'][:, plev, :].values, UM133ens_grid['lsrain'][n,:,0,:].values)
            #UM044_thetap = calc_thetap(UM044_plev['temp'][:, plev, :].values, UM044_plev['p'][plev].values,
            #                           UM044_plev['q'][:, plev, :].values, UM044ens_grid['lsrain'][n,:,0,:].values)
            P = str(UM133_plev['p'][plev].values)
            fname='Vid/thetap_plev_%03i.png'%tidx
        else:
            coldpool_co[i] = {}
            UM133_thetap = UM133_thetae[n,:,0]
            UM044_thetap = UM044_thetae[n,:,0]
            fname='Vid/thetap_surf_%03i.png'%tidx
            P = 'Surface'
        for j, (name, data, coords) in enumerate(zip(('UM133', 'UM044'), (UM133_thetap, UM044_thetap), (UM133_rain, UM044_rain) )):
            tit = ('UM 1.33km','UM 0.44km')[j]+ r' (\rom{{{}}})'.format(i)
            loc = comp[i][name]['mean'].values.argmax()
            #tidx = idx[i][name]['time'][min(loc,len(idx[i][name]['time'])-1)]
            
            time = pd.DatetimeIndex(UM044_plev.coords['t'].values).round('10min')
            T1 = pd.DatetimeIndex(UM044ens_grid.coords['t'].values).round('10min')[tidx]
            dt = np.fabs((time - T1).total_seconds()).argmin()
            locidx = idx[i][name]['loc'][loc]
            
            
        
            ilon=0,-1
            ilat=0,-1
        
            lons = coords.coords['longitude'][ilon[0]:ilon[1]].values
            lats = coords.coords['latitude'][ilat[0]:ilat[1]].values
            tsurf = data[dt, ilat[0]:ilat[1], ilon[0]:ilon[1]]
            tsurf = tsurf-tsurf.mean()
            minloc = np.unravel_index(tsurf.argmin(), tsurf.shape)
            maxloc = np.unravel_index(tsurf.argmax(), tsurf.shape)
            
            #grad = np.gradient(data[tidx, pilat[0]:pilat[1], pilon[0]:pilon[1]])
            #fulgrad = np.sqrt(grad[0]**2 + grad[1]**2)
            #peak = np.unravel_index(fulgrad.argmax(), fulgrad.shape)
            if first:
                pax = fig.add_subplot(2,2,nplot)
                m = Basemap(llcrnrlat=min(lats), llcrnrlon=min(lons), urcrnrlat=max(lats), urcrnrlon=max(lons),
                            resolution='f', area_thresh=1, ax=pax)
                im = m.pcolormesh(lons, lats, tsurf, vmin=-3, vmax=3, cmap=colmap.GMT_polar, shading='gouraund')
                M.append((pax, m, im))
            else:     
                pax, m, im = M[nplot-1]
                im.set_array(tsurf.ravel())
            m.drawcoastlines(linewidth=0.5)
            pax.tick_params(labelsize=12)
            
            
            if pdata == 'pres':
                bla = 'surf'
                #coldpool_co[i][name] = (111, (lons[32], lats[10]), (lons[23], lats[20]))
                if n == 0 and j == 1:
                    coldpool_co[i][name] = (111, (130.92701944, -11.721090399999999), (lons[maxloc[1]]-0.03, lats[maxloc[0]]+0.1))
                else:
                    coldpool_co[i][name] = (111, (130.92701944, -11.721090399999999), (lons[maxloc[1]], lats[maxloc[0]]))
                #coldpool_co[i][name] = (111, (lons[minloc[1]], lats[minloc[0]]), (lons[maxloc[1]], lats[maxloc[0]]))
                
                #(tidx, (locidx[0], locidx[1]), (speak[0], speak[1]) 
            tidx, locidx, speak = coldpool_co[i][name]

            m.scatter(locidx[0], locidx[1], s=20, c='r')
            m.scatter(speak[0], speak[1], s=20, c='g')
            pax.set_title(tit, fontsize=fontsize)
            nplot += 1
    if first:
        fig.subplots_adjust(right=0.98, bottom=0.25, top=0.99,left=0.01, hspace=0.2, wspace=0.15)
        cbar_ax = fig.add_axes([0.14, 0.20, 0.74, 0.01])
        cbar=fig.colorbar(im, cax=cbar_ax, orientation='horizontal')
        cbar.ax.tick_params(labelsize=fontsize-2)
        tlable = time.tz_localize(utc).tz_convert(timezone).tz_localize(None)[dt].strftime('%d/%m %H:%M')
        cbar.set_label('%s hPa Density Potential Temp. Pertubation at %s [$^\circ$C]' %(P,tlable),size=fontsize)
    #plt.suptitle('Time: %s LT'%time[0] .strftime('%Y-%m-%d %H:M'))\
    #fig.savefig(fname, bbox_inches='tight', format='png', dpi=72)
    first = False
    
dest_dir = os.path.abspath('Vid')
#fig.clf(), plt.close()
#sys.stdout.write('... ok\n')
#make_mp4_from_frames('Vid', dest_dir, 'Thetae_surf.mp4', 4, glob='Vid/thetae_*.png')
In [1429]:
from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes, inset_axes
from mpl_toolkits.axes_grid1.inset_locator import mark_inset
with xr.open_dataset(os.path.join(UMdir, 'tmp', 'ancil', 'qrparm.orog.nc')) as tmp:
    oro = tmp.variables['ht'][0,0]
    lats = tmp.variables['rlat'][:]
    lons = tmp.variables['rlon'][:]
    fig = plt.figure()
    ax = fig.add_subplot(111)
    m = Basemap(llcrnrlat=min(lats), llcrnrlon=min(lons), urcrnrlat=max(lats), urcrnrlon=max(lons),
                resolution='f', area_thresh=1, ax=ax)
    im = m.pcolormesh(lons, lats, oro, vmin=-50, vmax=110, cmap=colmap.GMT_globe, shading='gouraud')
    m.scatter(lons[161], lats[83], s=200, c='r')
    m.scatter(lons[141], lats[63], marker='*', s=100, c='k')
    m.scatter(lons[181], lats[103], marker='*', s=100, c='k')
    
    m.plot([lons[141], lons[181]], [lats[63], lats[103]])
    locidx = list(bresenham(141, 63, 181, 103))
    ht = get_data2(oro, locidx)
    m.drawcoastlines()
    #ax = fig.add_subplot(212)
    cd1 = dist((lons[161], lats[83]), (lons[181], lats[103]))
    cd2 = dist((lons[141], lats[63]), (lons[161], lats[83]))
    dd = list(-np.linspace(0,cd2,len(ht)//2)[::-1])+list(np.linspace(0,cd1, len(ht)//2 + 1))
    axins = inset_axes(ax, width='25%', height='25%', loc=2)
    axins.fill_between(dd, ht)
    axins.grid()
    axins.set_xlim(min(dd), max(dd))
    axins.set_title('Top at %0.2fE/%0.2fS'%(lons[161], np.fabs(lats[83])), fontsize=32)
In [1408]:
ax.fill_between?
  • Cross-section of that snapshot
In [1130]:
fig=plt.figure(figsize=(15,10))
fontsize=16
ds = 10
hfmt = dates.DateFormatter('%H:%M')

lim =800
for dt in (17,):
    nplot = 1
    for n, i in enumerate((3, 5)):
        UM133_plev, UM044_plev = get_file('vert_cent', num=i-1, UMdir=os.path.join(UMdir, '1p33km'),
                 first='20061112_0000', last='20061112_1200', remap_res='1.33km')
        UM133_rain, UM044_rain = get_file('rain', num=i-1, UMdir=os.path.join(UMdir, '1p33km'),
                 first='20061112_0000', last='20061112_1200', remap_res='1.33km')
        P = UM133_plev['p'].values
        UM133_thetap = calc_thetap(UM133_plev['temp'][1:, :, :], P.reshape(1,-1,1,1),
                                   UM133_plev['q'][1:, :, :].values, UM133_rain['lsrain'][:-1,:].values)
        UM044_thetap = calc_thetap(UM044_plev['temp'][1:, :, :], P.reshape(-1,1,1),
                                   UM044_plev['q'][1:, :, :].values, UM044_rain['lsrain'][:-1,:].values)
        inp = (UM133_thetap, UM044_thetap)
        lookup={1:UM044_plev, 0:UM133_plev}
        '''
        if n == 1:
            inp = inp[::-1]
            lookup={0:UM044_plev, 1:UM133_plev}
        else:
            lookup={0:UM133_plev, 1:UM044_plev}
        '''
        for j, (name, data) in enumerate(zip(('UM133', 'UM044'), inp)):
        
            pax = fig.add_subplot(2,2,nplot)
            tit = ('UM 1.33km','UM 0.44km')[j]+ r' (\rom{{{}}})'.format(i)
            loc = comp[i][name]['mean'].values.argmax()
        
            tidx2, co1, co2 = coldpool_co[i][name]
            tidx2 = 111
            co_data = lookup[j]
            try:
                lats = co_data.coords['latitude'].values
                lons = co_data.coords['longitude'].values
            except KeyError:
                lats = co_data.coords['lat'].values
                lons = co_data.coords['lon'].values
            time = pd.DatetimeIndex(co_data.coords['t'].values).round('10min')
            T1 = pd.DatetimeIndex(UM133ens_grid.coords['t'].values).round('10min')[tidx]
            gdist = dist(co1, co2)
            time = time.tz_localize(utc).tz_convert(timezone).tz_localize(None)[dt]
            #dt = np.fabs((time - T1).total_seconds()).argmin()
            co1 = np.fabs(lons-co1[0]).argmin(), np.fabs(lats-co1[1]).argmin()
            co2 = np.fabs(lons-co2[0]).argmin(), np.fabs(lats-co2[1]).argmin()
            locidx = list(bresenham(co1[0], co1[1], co2[0], co2[1]))
            while len(locidx) < 6:
                xinc = np.sign(locidx[-1][0] - locidx[0][0])
                yinc = np.sign(locidx[-1][1] - locidx[0][1])
                locidx = list(bresenham(co1[0], co1[1], locidx[-1.][0]+xinc, locidx[-1][1]+yinc))
            a= np.array(locidx)
       

            thetae = get_data2(data[dt] - 273.15, locidx)
            X= np.linspace(0,gdist,thetae.shape[0])
            p = P[P>=lim]
            xlim = len(p)
            thetae=thetae.T[:xlim]
            p1 = np.linspace(lim,1000, 200)[::-1]
            p2 = P[:xlim]
            x1 = X
            th1 = thetae
            ds = xr.Dataset({'thetae': (('p', 'x'), th1)}, {'x': X, 'p': p2})
            ds = ds.interpolate_na(xc=x1, yc=p2, method='polynomial', order='12')
            th1 = ds['thetae'] - ds['thetae'].mean()
            th1 = th1[:,::-1]#/th1.mean()
            im = pax.pcolormesh(ds.coords['x'], ds.coords['p'], th1, vmin=-18/9.81, vmax=18/9.81, cmap=colmap.GMT_haxby, shading='gouraud')
            #pax.scatter(X[peak[0]], P[1], c='k')
            if nplot in (1, 3):
                pax.set_ylabel('Pressure [hPa]', fontsize=fontsize)
            else:
                pax.axes.yaxis.set_ticklabels([])
            if nplot in (3, 4):
                pax.set_xlabel('Distance [km]', fontsize=fontsize)
            else:
                pax.axes.xaxis.set_ticklabels([])
            pax.invert_yaxis()
            pax.set_ylim(1000,lim)
            pax.set_xlim(0,40)
            pax.tick_params(labelsize=fontsize-2)
            pax.set_title(tit, fontsize=fontsize)
            nplot += 1

    fig.subplots_adjust(right=0.98, bottom=0.08, top=0.99,left=0.01, hspace=0.2, wspace=0)
    cbar_ax = fig.add_axes([0.14, 0.0, 0.74, 0.017])
    cbar=fig.colorbar(im, cax=cbar_ax, orientation='horizontal')
    cbar.ax.tick_params(labelsize=fontsize-2)
    cbar.set_label('Density Potential Temp. Pertubation at %s LT [$^\circ$C]'%(time.strftime('%H:%M')),
                   size=fontsize)
    #fig.savefig(os.path.join(outdir,'crossect_%03i.png'%dt), bbox_inches='tight', format='png', dpi=72)
#fig.clf()
#plt.close()
#make_mp4_from_frames(outdir, dest_dir, 'Corssect_nativ_2.mp4', 2, glob='crossect_*')
In [353]:
tmp_dir='/home/unimelb.edu.au/mbergemann/Data/Extremes/UM/darwin/RA1T/1p33km'
fluxes={3:{}, 5:{}}
for en in native.keys():
    tmp044, tmp133 = get_file('vert_cent', num=en-1, UMdir=tmp_dir, first='20061112_0000', last='20061112_1200', remap_res='1.33km')
    rain044, rain133 = get_file('vert_wind', num=en-1, UMdir=tmp_dir, first='20061112_0000', last='20061112_1200', remap_res='1.33km') 
    
    fluxes[en] = {'UM133': (tmp133, rain133), 'UM044': (tmp044, rain044)}
  • The day (ensemble member III and V) in their native resolutions:
In [1131]:
embed_vid(os.path.join('Vid','ColdPool_nativ_2.mp4'))
Out[1131]:
In [1132]:
embed_vid(os.path.join('Vid','WeekOfHector-Ens-2.mp4'))
Out[1132]:
  • Distributions of some variables
In [1377]:
mpld3.disable_notebook()
variables=dict(omega=('Vert. motion', (0, 10)), bflux=('Bouyancy Flux', (-3,30)), cloud_z=('Cloud Top', (4,12)),
               momflux=('Momentum Flux', (0., 0.3)), moistflux=('Moisture Flux', (-0.01, 0.04)))
               #cloud_w=('Cloud Water', (0,0.005)))
fig = plt.figure()
fig.subplots_adjust(right=0.94, bottom=0.025, top=0.6,left=0.01, hspace=0, wspace=0)
i = 0
for y, prop in variables.items():
    label, ylim = prop
    npl = len(list(variables.keys()))
    i += 1
    ax = fig.add_subplot(npl,1,i)
    ax = sns.boxplot(x="quant", y=y, hue="run", data=storm_prop, palette="muted", ax=ax)
    #ax = sns.stripplot(x="quant", y=y, hue="run", data=df_stack, jitter=True, palette="Set2", dodge=True)
    if i == 1:
        ax.legend(loc=2, fontsize=24)
    else:
        ax.legend(fontsize=1)
    ax.tick_params(labelsize=24)
    ax.set_xlim(0.5,5.5)
    ax.set_ylim(*ylim)
    ax.yaxis.set_ticks(ax.yaxis.get_ticklocs()[1:])
    ax.set_xlabel('Quintiles', fontsize=26)
    ax.set_ylabel(label, fontsize=26)
  • Some vertical profiles
In [1378]:
fig = plt.figure()
from matplotlib.ticker import Formatter

from scipy import interpolate
variables=dict(bflux=('Bouyancy Flux', (-0.2,1.02)), momflux=('Momentum Flux', (2e-4, 9e-2)),
               #cloud_z=('Cloud Water/Ice', (0,2.4e-4)),
               mflux=('Moisture Flux', (0, 8e-4)), omega=(('Vert. Motion'), (-9.5e-2, 4.5e-1)))
i = 0
fontsize=24
y = np.linspace(P[-1], P[0], 40)[::-1]
for x, prop in variables.items():
    label, xlim = prop
    npl = len(list(variables.keys()))
   
    for ii in range(1, 6):
        i += 1
        ax = fig.add_subplot(npl,5,i)
        for ff, fluxes in enumerate((fluxes044, fluxes133)):           
            flux = fluxes[x][ii]
            f = interpolate.interp1d(P, flux, kind='slinear')
            #xx = f(y)
            y = P
            xx = flux
            ax.plot(xx, y, label=('UM 0.44km', 'UM 1.33km')[ff])
            ax.set_ylim(max(y), min(y))
            ax.set_xlim(*xlim)
            if ii == 1:
                ax.set_ylabel('Pressure [hPa]', fontsize=fontsize-2)
                ax.legend(loc=1, fontsize=fontsize)
            else:
                ax.yaxis.set_ticklabels([])
            ax.set_xlabel(label, fontsize=fontsize-2)
            ax.tick_params(labelsize=fontsize-2)
            ax.ticklabel_format(axis='x', style='sci')
            if i < 6:
                ax.set_title('Qantile  %i'%ii, fontsize=fontsize)
            #ax.xaxis.set_major_formatter('%.1E')
            ax.ticklabel_format(style='sci', axis='x') 
            ax.xaxis.major.formatter.set_powerlimits((0,0))
fig.subplots_adjust(right=0.98, bottom=0.25, top=0.99, left=0.01, hspace=0.15, wspace=0.04)
In [1125]:
calc_thetav??